#!perl

use strict;
use Getopt::Long;
use MaterialsScript qw(:all);


# change the name of input and output file, and threshold
my $doc= Documents->Import("Layer1.xtd");
my $threshold = 26;
my $file = Documents->New("Hbonds.txt");


#------------------------------------------------
my $NumberOfMolecules = $doc->UnitCell->Molecules->Count;
my $trajectory = $doc->Trajectory;

for (my $frame=1; $frame<=$trajectory->NumFrames; ++$frame){
	$trajectory->CurrentFrame = $frame;
	my $newfile = Documents->New("temp_file.xsd");
	$newfile->CopyFrom($doc);
	foreach my $molecule (@{$newfile->UnitCell->Molecules}){
		if ($molecule->chemicalFormula ne "H2 O"){
			$molecule->Delete;
		}
		 
		if ($molecule->center->Z < $threshold){
			$molecule->Delete;
		}
	}
	$newfile->CalculateHBonds;
	my $HBonds = $newfile->UnitCell->HydrogenBonds;
	my $NumberofHBonds = $HBonds->Count;
	$file->Append(sprintf "%d,%d\n", $frame,$NumberofHBonds);
	$newfile->Delete;
}
